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^: ABSTRACT 

t> ■ The structure of the coronal magnetic field prior to eruptive processes and the 

QQ ■ conditions for the onset of eruption are important issues that can be addressed 

^ ■ through studying the magnetohydrodynamic stability and evolution of nonlinear 

/ ■ force-free field (NLFFF) models. This paper uses data-constrained NLFFF mod- 

O ■ els of a solar active region that erupted on 2010 April 8 as initial condition in 

MHD simulations. These models, constructed with the techniques of fiux rope 
insertion and magnetofrictional relaxation, include a stable, an approximately 
k^ . marginally stable, and an unstable configuration. The simulations confirm previ- 

;-H ■ ous related results of magnetofrictional relaxation runs, in particular that stable 

■ ■ ■ fiux rope equilibria represent key features of the observed pre-eruption coronal 

structure very well and that there is a limiting value of the axial fiux in the rope 
for the existence of stable NLFFF equilibria. The specific limiting value is lo- 
cated within a tighter range, due to the sharper discrimination between stability 
and instability by the MHD description. The MHD treatment of the eruptive 
configuration yields very good agreement with a number of observed features like 
the strongly inclined initial rise path and the close temporal association between 
the coronal mass ejection and the onset of fiare reconnection. Minor diflFerences 
occur in the velocity of fiare ribbon expansion and in the further evolution of 
the inclination; these can be eliminated through refined simulations. We suggest 



that the shngshot effect of horizontaUy bent flux in the source region of eruptions 
can contribute signiflcantly to the inchnation of the rise direction. FinaUy, we 
demonstrate that the onset criterion formulated in terms of a threshold value 
for the axial flux in the rope corresponds very well to the threshold of the torus 
instability in the considered active region. 

Subject headings: magnetohydrodynamics (MHD) — Sun: corona — Sun: coro- 
nal mass ejections (CMEs) — Sun: fllaments — Sun: flares — Sun: magnetic 
flelds 



1. Introduction 



The coronal source regions of solar eruptions (eruptive prominences, coronal mass ejec- 
tions, and flares) are characterized by low plasma beta, especially in active regions, where 
/3 ~ 10~^ (lGarvll200ll ). Consequently, the magnetic fleld must remain nearly force free in the 
equilibrium sequence seen as the quasi-static evolution toward the onset of eruption . Since 
the e nergy needed to power the eruption can only be stored in coronal currents ( JForbes 
2000) and s ince these curre nts tend to be strongly concentrated in the core of the source 
region (e.g. JSun et al.ll2012l ). the pre-eruptive fleld is, to a good approximation, a nonlinear 
force-free fleld (NLFFF) obeying V X S = cxB with the "force-free parameter" a{x) being 
a scalar function that varies across B. 

In order to understand the mechanism of eruptions and to quantify the criteria for 
their onset, knowledge of the fleld in the erupting coronal volume is required. However, the 
three-dimensional (3D) fleld structure is not amenable to measurement, and reliably infer- 
ring the coronal NLFFF from photospheric magnetograms has proven very diflicult, even 
if vector magnetograms are available. Two strategies have been pursued. Extrapolation 
techniques solve the force-free equation numerically using a vector magnetogram as bottom 
boundary condition. Several schemes produced reliable results when applied to test flelds, in 
particular but not exclusivel y, the Grad-Rubin ite ration, an optimization scheme, and mag- 
netofrictional relaxation; see lSchrijver et al.l ( 120061 ) for an overview. However, it appears that 
the reconstruction of the coronal fleld from observed vector magne tograms matches realit y 



closely only in a s till re la tively small number of c ases; se e especially ISchriiver et a" 



Canou fc Amaril fcoioh . lYelles Chaouche et al.l f l2012h . IValori et al.l ( 120121 ). and 



(l2008ah. 



Sun et al. 



(l2012r). The methods had diflSculty in producing reliable resul ts for some other cases (e.g 



Metcalf et al 



2008 



Schrijver et al.ll2008al : iDe Rosa et al.ll2009l ). The violation of the force- 



free condition at the photospheric level in part of the magnetogram, uncertainties of the 
transverse magnetogram components, and the fragmentation of the flux in the photosphere 



down to sub-resolution scales are supposed to cause the problems which are not yet clearly 
understood. 



The ex trapolation technique 
instability (iXorok fc Klieml boOS : 



las supported the mo deling of eruptions as a flux rope 



Kliem fc Torokl 120061 ) by computing stable and unsta- 



ble coronal fields ( IValori et al.l l2010l ) from magnetograms of the active-region model by 
Titov fc DemoulinI ( 119991 ). However, to d ate only few extr a polations of obse r ved m agne- 
tograms have found a flux rope, see, e.g., lYan et al.l ( 120011 ). ICanou fc Amaril (120101 ). and 
Yelles Chaouche et al.l ( 120121 ) . A counterexample is the reconstruction of a highly nonpoten- 



tial active region on the Sun (JNindos et al.l l2012l : Valori, private comr aunication) , although 



the subsequent viole nt eruption was suggestive of a flux rope eruption ( ISchrijver et al. 
Zharkovet aPboilh . 



2011 



The flux rope insertion method (Ivan BallegooijenI I2004J : Ivan Ballegooijen et al.l l2007l ) 
provides a viable alternative tool for the determination of the coronal NLFFF. Here a flux 
rope is inserted in the potential fleld computed from the vertical magnetogram component, 
to run at low height along the magnetogram's polarity inversion line (PIL) in the section 
where a fllament channel indicates nearly horizontal, cur rent-carrying fleld . This conflgu- 
ration is then numerically relaxed using magnetofriction ( JYang et al.lll986l ). The resulting 
numerical equilibrium has received observational support in a growing number of cases by 
favorable comparison of various fleld line shapes (dips, arches of various shea r angle rela- 
tive to the PIL, S shape) and quasi-separatrix layers with coronal structures (iBobra et al. 



20081: 



Su et al. 



Savcheva et al. 



20091. 



2012c 



201ll : ISu fc van BallegooijenI l2012l : ISavcheva fc van BallegooijenI l2009l : 
. In these studies, quiet-Sun or decaying active regions were mod- 



eled, where persistent shearing and flux cancelation make the formation of a coronal flux 
rope likely. It should be noted, however, that the method is not restricted to conflgurations 
containing a fl ux rope; the numerical relaxation can also result in an arcade conflguration 
( ISu et al.ll201ll ). The application of the flux rope insertion method led to the suggestion that 
the ratio between the magnetic flux in the rope, primarily the axial flu x, and the unsigne d 
flux in the active region may provide a new onset criterion for eruptions ( JBobra et al.ll2008l) 
This has found considerable support in subsequent appli cations (see in particular ISu et al. 



2011) as well as in MHD simulations of flux cancelation (lAulanier et al. 



2010l ). although the threshold may vary in a wide range ( ISavcheva et al. 



2010l: Amari et al. 



2012bh . While the 



magnetofrictional relaxation method is well suited to flnd the dividing line between stabil- 
ity and instability in parameter space, it fails in correctly describing the dynamic evolution 
of unstable conflgurations because it represents a strongly reduced MHD description which 
excludes the MHD waves. It is therefore of interest to study the relaxed or partially re- 
laxed conflgurations in full MHD simulations. This allows further judgment of the obtained 
conflgurations by comparison with observations of erupting regions, and, additionally, an 



independent test of the marginal stability line in parameter space. The present investigation 
aims at carrying out the first such experiments. 



Su et al.l ( 1201 ll . hereafter Paper I) performed extensive and detailed modeling of NOAA 



Active Region (AR) 11060 at the stage just prior to its eruption on 2010 April 8. The event 
included an erupting filament that evolved into a coronal mass ejection (CME) accompanied 
by a moderate flare. By varying the parameters of the inserted flux rope in a range of axial 
($axi) and poloidal (Fpoi) flux, the stability boundary in the $axi~^poi plane was determined 
with relatively high accuracy. Weakly twisted and arching field lines of a resulting stable 
configuration rather close to the stability boundary showed good correspondence to X-ray 
and EUV loops prior to the eruption. Low-lying highly sheared field lines of an unstable 
configuration rather close to the boundary corresponded well to the observed initial flare 
loops. A further result of high interest with regard to the onset of eruptions was that an 
X-type magnetic structure known as a hyperbolic flux tube (JTitov et al.l l2002l : iTitovl 120071 ) 
appeared at low heights under the flux rope when the axial flux was raised to approach the 
stability boundary. In the following we study a series of configurations from this investigation 
in the zero-beta MHD approximation. The configurations are arranged on a path in the 
$axi-^poi plane that crosses the stability boundary and include the models which compare 
favorably to the observations. 

The simulations also allow us to compare the condi tions at the onset of instability with 
the ratio of rope an d ambient flux (IBobra et al.l l2008r) and the heig ht profile of the flux 
rope's ambient field (Ivan Tend fc Kuperuslll978l : lKliem fc Torokl 12009 ). Identification of the 
conditions necessary for instability is a key goal of this research. 



Section [2] presents a summary of the observations from lSu et al.l (J201ll ). This is followed 
by a detailed description of the numerical model in Sections [3] and [H The relaxation of stable 
configurations and the eruption of an unstable configuration are presented in Sections [5] and 
[6l respectively. Our conclusions are given in Section [71 



2. Summary of Observations 

The event under study originated in AR 11060 on 2010 April 8 around 02:30 UT. 
It included an eruptive prominence that evolved into a CME, a two-ribbon solar flare of 
GOES class B3.7, prominent coronal dimmings, and a coronal E UV wave. The erup tion was 
observed on disk by the Atmospheric Imaging Assembly (AIA, iLemen et al.ll2012r) on the 
Solar Dynamics Observatory (SDO) and by the X-Ray Telescope (XRT, iGolub et al.ll2007l ) 
on the Hinode mission. The Extreme Ultraviolet Imager (EUVI) and the COR coronagraphs 
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Fig. 1. — Structure of AR 11060 before and during the eruption observed by (a) Hin- 
ode/XRT^ and (b~d) SDO/AIA. In panel (a) the contours refer to the photospheric hne of 
sight (LOS) magnetogram taken by SDO/HMI at 02:00 UT on 2010 April 8 (white: positive; 
black: negative). The coronal loops in the core of the region are organized in three sets, 
labeled as 1-3 and discussed in the text. (A color version of this figure is available in the 
online journal.) 



(JHoward et al.ll2008l ) aboard the Solar-Terrestrial Relations Observatory (STEREO) Ahead 
spacecraft saw the event nearly exactly at the east limb. A detailed description of these data 
including animatio ns can be found in Paper I. Additional information is given in the analysis 
of the EUV wave (JLiu et al.l l2010l ) and of the possible occurrence of the K elvin-Helmholtz 
instability at the surface of the expanding flux (jOfman fc Thompsonll201ll ). 



Figure [I shows Hinode/XRT (a) and SDO/AIA (b) images of AR 11060 before the 
eruption. The core of this active region contains three sets of coronal loops as shown in 
Figure [11(a). The set of loops labeled 1 is located in the northwestern part of the region 
(white arrow). Loop sets 2 and 3 are associated, respectively, with the major positive 
polarity (white contours) in the northeast of the region and with a group of minor positive 
flux patches in the southeast (both marked by black arrows). The combined loop sets 1 



and 2 yield a slightly sigmoidal appearance at soft X-rays. Figure [T](b) shows a thin, dark 
filament following mainly the highly sheared loop set 2, i.e., located between loop sets 1 and 
3. 

The eruption began around 02:10 UT with motion of material along this filament in 
the southeastward direction (as seen by AIA) and nearly horizontally (as seen by EUVI-A 
on STEREO Ahead). About 18 minutes later, the filament began to lift off. From about 
02:33 UT, the overlying loops are seen to rise in the EUVI-A 195 A images. The first fiare 
brightenings appear in AIA 94 A images on either side of the erupting filament around 
02:30 UT, and the GOES soft X-ray emissions commence near 02:45 UT. Thus, the onset of 
upward fiux expansion which evolved into the CME and the onset of the fiare brightenings 
occurred in close temporal and spatial association. 

Initially highly sheared fiare loops became visible between 02:30 and 02:40 UT in the 
AIA 94 A channel (Figure [T](c)), they connected the first fiare brightenings in the strong-field 
section of the PIL. The typical evolution toward lower shear was seen in the loops forming 
subsequently (Figure [T](d)). These observations indicate that the component of the magnetic 
field along the PIL points northwestward and that the field has right-handed helicity, consis- 
tent with the forward S shape of the combined loop sets 1 and 2. The brightenings early in 
the fiare (Figure [T](c)) as well as the subsequent equatorward growth of the CME dimmings 
yield an overall inverted U shape of the parts of the region actively involved in the eruption. 
This is similar to the combined loop sets 1 and 3 and suggests that their fiux is an integral 
part of the eruption, although loop set 3 connects to a minor fiux area in the magnetogram. 

The erupting fiux took a strongly inclined path toward the equator, initially at about 
45° from vertical, then developing an even more inclined southward expansion in the low 
corona imaged by EUVI-A, and finally turning into a more radial propoagation near the 
equatorial plane, as imaged by the C0R2-A coronagraph. The EUVI data do not reveal the 
rise profile of the erupting filament, as the material very quickly lost contrast in the 195 A 
channel and the cadence of the other channels was too low. However, the rise velocity of 
overlying coronal loops can be estimated and is found to reach about 170 km s~^ within 
the EUVI-A field of view (see Section [6Tl) . The CME reached a median projected speed 
of ^ 520 km s-Mn th e C0R2-A field of view, as quoted in the CACTus CME Catalog 



(IRobbrecht et al.ll2009h . 



^http://secchi. nrl.navy.mil/cactus/ 



3. Magnetofrictional Modeling 

The presence of highly sheared loops (Figured]) indicates that the coronal magnetic field 
in the observed region deviates significantly from a current-free potential field. In Paper I we 
used a set of codes, called the Coronal Modeling System (CMS), to construct non-potential 
field models of the observed region. The methodology involves inserting a thin magnetic fiux 
rope into a three-dimensional (3D) potential -field model along a specified path, and then 
applying magnetofrictional relaxation (MFR) ( JYang et al. II 19861 : Ivan Ballegooijen et al.ll2000l ) 



to produce either a non-linear force-free field (NLFFF) model or unstable model, depending 
on the magnitude of the inserted fiux. MFR refers to an evolution of the magnetic field 
B{r^t) in which the medium is assumed to be highly conducting and the plasma velocity 
V is proportional to the Lorentz force j X B, where j = V X S is the current density. 
The resulting 3D magnetic models are based on observed photospheric magnetograms and 
therefore accurately represent the lower boundary conditions on the magnetic field in the 
fiaring active region. This makes such models ideally suited as initial conditions for 3D MHD 
simulations. 

In this paper we use fiux rope insertion and MFR to produce initial conditions for a 3D 
MHD code, which is described in Section [H One problem is that the MHD code assumes 
Cartesian geometry, whereas the CMS codes use spherical geometry. Therefore, the methods 
used in Paper I need to be modified in several ways, as described in the following. First, 
we construct a map of Br{9^ 0), the radial component of magnetic field at the solar surface 
(r = Rq) in and around the observed active region. The map is based on an SDO/HMI 
magnetogram taken at 2:00 UT. The map is computed by interpolating the observed LOS 
field onto a longitude-latitude grid, and then dividing by the cosine of the heliocentric angle 
to obtain B^. This original map has 384 x 384 cells and covers an area of 41.3° in longitude 
by 36.8° in latitude, centered on the observed active region. Then the spatial scale of the 
map is reduced, and the center of the map is displaced to the equator, corresponding to a 
reduction in scale by a factor / = 16.97. The rescaled map covers only 2.2° in longitude and 
latitude, which is sufficiently small that the resulting spherical models can be treated as if 
they were Cartesian. 

In CMS the magnetic field is described using two spatial grids: a high resolution grid 
(HIRES) covering the target region and its local surroundings, and a low resolution grid 
covering the entire Sun. The non-potential field in the HIRES region is described in terms 
of vector potentials (S = V X A). In the standard approach, the global grid is used for a 
potential field source surface (PFSS) model; its main purpose is to improve the side boundary 
conditions for the HIRES domain. However, for the scaled models used here the global grid 
can no longer be used. Instead we assume that the normal component of magnetic field 




Fig. 2. — Ha image of the observed active region, obtained at Kanzelhohe Solar Observatory 
on 2010 April 8 at 9:00 UT. The red and green contours show the photospheric magnetic 
field as derived from the SDO/HMI magnetogram taken at 2:00 UT. The blue curve shows 
the path of the inserted flux rope between the marked end points. (A color version of this 
figure is available in the online journal.) 

vanishes at the side boundaries of the HIRES domain {Bo = and 5^ = at the latitudinal 
and longitudinal boundaries, respectively). These boundary conditions imply that the net 
flux entering the domain through the lower boundary is also present at larger heights, and 
leaves the domain through the outer boundary. We found that small flux imbalances between 
positive and negative magnetic fluxes in the imposed magnetic map {B^) can produce flelds 
at large heights that are dominated by such "monopole" components. Furthermore, in 
preliminary MHD simulations we found that such monopolar flelds can inhibit the eruption 
of low-lying flux ropes. Therefore, in the present paper we correct the flux imbalance in 
Br{9^(j)) at the lower boundary. The photospheric fluxes are balanced by subtracting a 
spatially constant value of about 2.2 Gauss from 5^(0,0), except near the border of the 
computational domain where S^ is set to zero. 

The next steps are to select the path of the flux rope, compute a potential fleld, and 
insert the flux rope into the 3D magnetic model. The path should follow the polarity inversion 
line of the active region; we use the same path as in Paper I. Figure [2] shows the selected 
path superposed on a Ha image obtained at Kanzelhohe Solar Observatory at 9:00 UT. The 
red and green contours indicate the photospheric flux distribution Br{6^ (j)) before scaling. 

In this paper we construct three different models with the same flux rope path, height. 



and cross section, but with different values of the axial flux of the flux rope. The axial fluxes 
before scaling are $axi = 4 x 10^°, 5 x 10^° and 6 x 10^° Mx, which were chosen because 
these values appear to straddle the stability boundary (see Paper I). After scaling these axial 
fluxes are reduced by the square of the scaling factor. The axial flux is inserted into the 3D 
model as a thin flux tube that is elevated above the photosphere. At the ends of the tube 
(indicated by circles in Figure [2j) , the flux of the tube is connected to surrounding sources 
on the photosphere. Poloidal flelds are added by wrapping flux rings around this tube. The 
assumed poloidal flux is Fpoi = 10^° Mx cm~^ (before reduction by th e scaling factor) . For 



more details on how the flux rope is inserted into the 3D models, see ISu et al.l (120111 ) and 
references therein. 

The next step is to apply 30,000 iterations of MFR to each of the three models. This 
drives the magnetic fleld toward a force- free state (if one exists) or causes the fleld to slowly 
expand (if the conflguration is unstable). Each 3D magnetic model is then resampled on 
the grid used by the MHD code. This grid is nonuniform in longitude, latitude and height; 
it provides high spatial resolution in the center of the active region and gradually reduced 
resolution farther away from the center. Figure [3] shows vertical cross-sections through the 
three models after resampling on the Cartesian grid. For consistency with the designation in 
Paper I we will refer to the axial flux values of the inserted flux ropes corresponding to the 
magnetogram size before the spatial scaling, 4 x 10^^, 5 x 10^^ and 6 x 10^^ Mx. The flrst 
panel indicates on a magnetic map at height 4.2 Mm where the vertical cross-sections are 
taken. The quantity plotted in the three other panels is a = j • S/S^, which is a measure of 
the parallel electric current. Note that the currents are concentrated near the edge of the flux 
rope, and that the height of the flux rope increases with axial flux. Higher axial flux implies 
higher total current running along the flux rope, which in turn i mplies a stronger repellin g 



force between the flux rope current and its subphotopheric image (JKuperus fc Raadulll974l ). 
so that the equilibrium height increases. Both models with $axi > 5 x 10^^ Mx possess an 
X-type magnetic structure — a hyperbolic flux tube — running at low height under the flux 
rope (see Paper I for detail). 



By applying the MFR much longer, ISu et al.l concluded in Paper I that, for the given 



value of Fpoi = 10^° Mx cm~^, the model with $axi = 5 x 10^^ Mx lies near the stability 
boundary in the $axi-^poi plane. All conflgurations with $axi ^ 4.5 x 10^° Mx relaxed to 
an apparently stable and approximately force-free state in 30,000 MFR iterations, while all 
conflgurations with $axi > 6 x 10^^ Mx did not relax, rather the inserted flux rope continued 
to rise. 
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Fig. 3. — The three configurations investigated in this paper after 30 000 MFR iteration 
steps. Shown are: (a) the magnetogram at a height of z = 4.2 Mm, the PIL at this height, 
and the position of a vertical plane used in the analysis; the plane is placed perpendicular 
to the PIL such that it passes very near the apex point of the rope's magnetic axis for 
all configurations, (b)-(d) Force-free parameter a and in-plane magnetic field vectors in 
this diagnostics plane (for z > 4.2 Mm; s denotes horizontal distance). The configurations 
possess axial fiux of 4, 5, and 6 x 10^° Mx for (b)-(d), respectively, and identical poloidal 
fiux of 10^° Mx cm~^. Lengths are given in i?©//, i.e., they would correspond to a unit 
of IRq if scaled back to the original magnetogram size. Crosses show the initial position 
of the fiuid elements used as start points of the field line tracing in the subsequent plots 
and animations; the trajectories of these fiuid elements are followed through the simulations. 
Rainbow-colored crosses are placed inside the FR near its edge, and green crosses are placed 
in the inner part of the overlying fiux. The plus signs mark the apex point of the fiux 
rope's magnetic axis. The motion of this fiuid element is used to generate the rise profiles 
in Figures El [HI andfT4l (A color version of this figure is available in the online journal.) 
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4. Magnet ohydro dynamic Modeling 

In the MHD modeling, the configurations obtained after 30,000 magnetofrictional itera- 
tions are used as the initial condition, So(r), for the numerical integration of the compressible 
ideal MHD equations in the zero-beta limit. We set the initial velocity to zero, Uq = 0, and 
choose po{r) = 5o(r)^/^ as model for the initial density. This implies a gradual decrease of 
the initial Alfven velocity Va = 5o(2/xoPo)~''^^^ with distance from the strong field in the core 
of the active reg ion, such that theheight profile Va{z) approximately follows the empirical 
height profile in IVrsnak et al.l ( J2002L their Figure 5). Such a choice of po{r) did generally 
yield good quantitative niatches with observed CME rise profiles in previous simulations 
(e.g., iTorok fc Kliemlboosl : ISchrijver et"aDl2008bl : Ixiiem et al.ll2012h . 



The equations used and the numerical scheme are detailed in iTorok fc Klieml (120031 ). 
The condition (3 = decouples the energy equation from the system. Gravity can be 
neglected for both purposes of the simulations in this paper: further relaxation of stable 
force-free configurations and the study of CME acceleration in a region of relatively strong 
magnetic field (an active region) in the case of instability. Numerical diffusion breaks the 
frozen-in condition when thin layers of inhomogeneous field (i.e., high current densities) 
develop, thus allowing magnetic reconnection. The equation of motion includes a viscous 
term solely to ensure numerical stability. 

Careful attention is paid to the amount of numerical diflFusion in the integration. We 
use a modified verson of the two-step Lax-Wendroff scheme that replaces the stabilizing but 
very diflPusive Lax term in the auxiliary step o f the iteration by so-called artificial smoothing 
(ISato fc Hayashil Il979l : iTorok fc Klieml l2003l ) . While the Lax term replaces the value of 
the integration variable by the average of the values at the six neighboring grid points, 
the artificial smoothing replaces only a fraction of the value, a <C 1, by the average at 
the neighboring points. For the magnetic field, this averaging is switched oflF completely 
in the relaxation runs {as = 0) to minimize any slow diffusive change of the force-free 
equilibrium, improving the convergence toward the equilibrium. The eruption of the unstable 
configuration can be followed with the same setting; however, a parametric study of the field 
smoothing yields the most reliable numerical results and best numerical stability if a small 
level of smoothing, as ^ 10~^, is adopted in the volume under the rising fiux rope, where 
a vertical current sheet develops (see Section [6] for detail). This operation is not applied 
in the bottom layer of the box, {z = 0}, so that there is no diffusion of the field in the 
magnetogram plane. 

Minimizing the numerical diffusion in the momentum equation allows us to capture any 
residual forces that result from our initial conditions being out of equilibrium. We have 
experimented with the values of au and u (the coefficient of viscosity) , and chosen them as 
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close to the limit of numerical stability as reasonably possible, i.e., a lowering of either of 
them by a factor of 2 would lead to numerical instability in the course of the simulation run. 
Best numerical stability is obtained by smoothing the density at the same level. In the main 
volume of the box dp = (Ju = 0.005. A gradual enhancement of the smoothing toward the 
bottom boundary, by a factor two, is necessary to allow stable integration in the presence of 
the strong flux gradients in the magnetogram and their associated currents. Uniform small 
viscosity is chosen: u = 0.002 (after normalization). 

The integration volume corresponds to 0.72Rq x 0.72Rq x IRq when the magnetogram 
is scaled back by the factor / to the original edge length. This volume is discretized by a 
Cartesian grid with uniform resolution of O.OO2i?0 (1.9 arcsec) in the active region and its 
immediate surroundings (|x| < O.li?©, \y\ < O.li?©, z < 0.2Rq) and a gradually increasing 
spacing in the outer parts, reaching O.OO65i?0 at the side boundaries and O.O15i?0 at the top. 
Closed boundaries are implemented at the sides and top of the volume by keeping the velocity 
at zero, including the first inner grid layer, so that the field vector in the boundary layers 
does not change. At the bottom boundary the velocity is kept at zero in the magnetogram 
plane, ?x(x,y, 0,t) = 0, which ensures that the normal magnetogram component B^ keeps 
its initial values, but allows the transverse components to change, either to approach an 
NLFFF, or to evolve consistently with an eruption. 

The box height (li?©) is chosen as the length unit. The field strength is normalized by 
the peak value, max(5o) = 955 Gauss, which lies in the magnetogram plane. Velocities are 
normalized by the peak value of the initial Alfven velocity, Vao, whose location coincides 
with the location of the highest field strength. The initial Alfven velocity in the body of 
the fiux rope is about half this value. Time is normalized by the corresponding Alfven time 
ta = Rq/Vaq. 



5. Relaxation of Stable and Nearly Marginally Stable Configurations 

In our MHD description both models with $axi < 5 x 10^° Mx are found to relax to 
a stable force- free equilibrium (NLFFF). As expected, the relaxation proceeds faster and 
deeper for the model with $axi = 4 x 10^^ Mx; nevertheless, the initial configuration clearly 
includes a small level of residual forces and the fiux rope experiences further reconnection 
with the ambient field in the course of the relaxation. The overall evolution is characterized 
in Figure [H From the field line plot at t = it is clear that the eastern part of the inserted 
fiux rope has evolved significantly in the course of the MFR: the fiux connected from the 
end point of the path shown in Figure [2] to the main positive polarity in the northeast of 
the region (i.e., in the direction of loop set 2 in Figured]) now bulges out to the south in 
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Fig. 4. — Stages in the relaxation of the model with $axi = 4 x 10^^ Mx. Shown are field lines 
in vertical view and the magnetogram ^^(x, y) with the gray scale saturated at one tenth of 
the peak value. Rainbow-colored field lines lie in a fiux shell slightly inside the boundary of 
the fiux rope and green field lines visualize the inner part of the overlying fiux (see Fig [3]) . 
(A color version of this figure is available in the online journal.) 



the direction of loop set 3. Close similarity to loop set 1 on the western side was achieved 
with relatively little change. The strong evolution of the eastern part continues in the course 
of the MHD relaxation, driven primarily by reconnection with the minor positive polarity 
in the southeast of the region. This wraps fiux from the minor polarity around the fiux 
rope, passing northward under the rope and from there arching over the rope toward the 
main negative polarity in the west. The bulging of the southern part of the rope is thus 
considerably enhanced. As a consequence, the fiux rope reconnects with the perturbing fiux 
and its whole eastern elbow splits. 

This interaction with the ambient fiux rooted in the minor positive polarity is shown 
in more detail in Figure [5] and its accompanying animation. Further reconnection in this 
area subsequently simplifies the structure of the inserted fiux which returns to an essentially 
unsplit fiux rope, shorter on the eastern side, where it now displays closer similarity in shape 
with loop set 2. New low- lying connections run from the minor positive polarity to the 
eastern part of the main negative polarity under the main body of the relaxed fiux rope; 
these are similar to loop set 3. Note that these connections differ from the ones at t = 0. 
Their existence is not even indicated in the potential field. Additionally, new high-arching 
field lines extend from the minor positive polarity with a dominant east-west orientation. 
They have absorbed some of the twist of the inserted fiux. One can speculate that such 
higher field lines confine plasma of lower density, so that these connections should be less 
visible in X-ray and EUV images. Some irregularly shaped emission does exist in this area 
at the temperature (T ^ 1.5 MK) sampled by the AIA 193 A channel (Figure [l](b)). It 
consists of moss essentially cospatial with the minor fiux concentration and of weak diffuse 
emission slightly southward, which has a dominantly east- west orientation reminiscent of the 
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Fig. 5. — Same as Fig. [4] but with additional field lines traced from a set of fixed points 
located at z = in the minor positive flux concentration in the southeast of the active 
region. The strong interaction of this flux with the inserted rope leads to the southward 
bulging of the rope and to multiple reconnections involving the splitting and the subsequent 
recovery of a single rope. The upper row shows the cube above the selected magnetogram 
area in perspective views oriented in —x direction. (A color version and an animation of this 
figure are available in the online journal.) 



high-arching field lines seen in the simulation from t ^ 15ta onward. 

More quantitative information about the relaxation, displayed in Figure [6l is obtained 
by monitoring the position and velocity of the fluid element initially at the apex point of 
the rope's magnetic axis (marked by a plus sign in Figure [3j). We first note a minor but 
quick shift in apex height at the very beginning of the run (t < Sta), which potentially 
results from two effects. One reason lies in the difference between the codes in obtaining 
the current density. While the MFR code runs on a staggered grid, the MHD code iterates 
all variables at the same positions, leading to differences in derived quantities between the 
codes which increase with degrading spatial resolution. Since all models include a relatively 
thin layer of enhanced current density at the edge of the flux rope which is resolved by only 
^ 5-10 grid cells (Figure [3j), the Lorentz force differs somewhat between the MFR and MHD 
descriptions in this critical layer, and so does the equlibrium height of the flux rope. It is 
clear, however, that the overall difference in the volume of the flux rope must be moderate 
at the resolution employed: the height of the flux rope's apex increases by only 15% (from 
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Fig. 6. — Height-time and velocity-time profile showing the relaxation of the model with 
$axi = 4 X 10^^ Mx. The fluid element initially located at the estimated apex point of the 
flux rope's magnetic axis, shown in Fig. [3] by a plus sign, is followed. The solid lines show 
the distance d(t) from the initial projection of the apex onto the bottom plane and the 3D 
velocity of the fluid element, |t^(t)|. The dashed lines show the height above the bottom 
plane and u^ii). Dotted lines in the logarithmic plot show downward velocities (|2i^(t)| for 
u^ < 0). 



ho = 0.025 to h{t = 2) = 0.029). Stronger diflFerences exist in the flrst few grid layers above 
the magnetogram, where the strong fragmentation of the flux is not fully resolved, resulting 
in small volumes of high current density. Although the associated Lorentz forces do not have 
a coherent direction, they may still contribute to the change in equilibrium height. This 
contribution can be estimated through a comparison with the relaxation behavior of the 
potential fleld, calculated on the spherical grid of CMS and resampled on the Cartesian grid 
of the MHD code in the same manner as the three active-region models. The fluid element at 
the same initial position as the one considered in Figure [6] rises by only 2.1% to h = 0.0255, 
also experiencing most of this displacement (80%) within the flrst two Alfven times. Thus, 
the contribution from discretization errors in the current density near the bottom of the 
box appears to be only minor. Second, part of the initial rise may be due to the status of 
relaxation after the 30,000 MFR iterations. Although already deep, as the subsequent rapid 
decrease of the velocities shows, it is obviously not yet complete (see above). The low level 
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Fig. 7. — Stages in the relaxation of the model with $axi = 5 x 10^^ Mx in the same format 
as Figure [H (A color version and an animation of this figure are available in the online 
journal.) 



of viscosity and velocity smoothing in the MHD relaxation allows the small residual forces 
to build up a noticeable velocity in the motion to the stable equilibrium height. 

Subsequently, a further rise to h = 0.033 occurs which is associated with the reconnec- 
tions discussed above. The velocity of the fiuid element shows a nearly monotonic decrease, 
with a weak transient enhancement during t ^ (5-20)rA related to the southward bulging 
of the flux rope. Following the initial jump (i.e., after t ^ 2rA), the vertical velocity falls by 
more than three orders of magnitude to u^ ^ 10~^ Vao, which is a clear signature that the 
flux rope relaxes deeply. The magnetic energy in the box decreases by 2.7% in the course of 
the relaxation. 

The model with $axi = 5 x 10^^ Mx follows a rather similar evolution in the course 
of its relaxation, albeit on a much longer time scale and relaxing somewhat less deeply 
(Figures [71 and [8]). In particular, the reconnection of the inserted flux rope with the minor 
positive flux concentration in the southeast of the active region yields a similar bulging, 
splitting, and subsequent trend of recovery. The resulting relaxed state is largely similar, 
but here the flux rope remains partly split. Again, the major part of the relaxed rope 
matches the observed loop sets 1 and 2 relatively well. The split-off part of the rope and 
some of the new connections between the minor positive polarity and the main negative 
polarity correspond well to loop set 3. Finally, diffuse, high-arching loops, with a dominant 
east-west orientation, are formed south of the minor positive polarity, as in the model with 
less axial flux. A noteworthy difference is the higher initial displacement of the flux rope 
by ^ 50% from ho = 0.036 to /i(2) = 0.053 in the flrst ^ 2ta- It can be seen in Figure [3] 
that the resolution of the current layer at the surface of the flux rope in the two models is 
rather similar, so the stronger initial displacement must primarily be due to a less complete 
relaxation in the course of the MFR. 
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Fig. 8. — Height-time and velocity-time profile showing the relaxation of the model with 
$axi = 5 X 10^^ Mx in the same format as Figure [6l 



After about IOOta the residual velocities lead to a very slow rise of the flux rope position, 
hardly visible in the fleld line plots but obvious in Figure El It is diflicult to judge whether 
this evolution is part of the relaxation or due to of a diffusive drift of the equilibrium in 
the long computation comprising slightly over 10^ iteration cycles. Animated fleld line plots 
(Figure [71) show that the rope continues to evolve over this long time period (some of the fleld 
lines change their photospheric connections considerably), and on average the velocity of the 
rope continues to decrease (apart from episodic enhancements); both properties suggest that 
the slow drift may be part of the relaxation. Nevertheless, it is unclear whether continuing 
this run can lead to further understanding of this model free from numerical effects, so 
the run was terminated. Adopting the scaling of the simulation in Section 16. 1[ the nearly 
400 Alfven times of this relaxation correspond to up to 2.3 days in reality, longer than the 
equilibrium can be assumed to be static. By the end of the relaxation, the magnetic energy 
in the box has dropped by 4.5%. The model appears stable in the MHD treatment, but is 
doubtlessly much closer to the point of marginal stability than the previous model. 

To further test for the stability of this model, we have perturbed it at the times of 
apparently deepest relaxation, t = 123rA and t = 395rA. The flux rope was raised to 
h ^ 0.105 by prescribing upward velocities at the apex in a sphere of radius equal to the 
minor flux rope radius for 15rA. Slightly below this height the unstable third model starts its 
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Fig. 9. — Evolution of the model with $axi = 5 x 10^° Mx following a strong perturbation 
applied at a stage of deep relaxation {t = 123rA). Position and velocity of the fluid element 
at the apex of the flux rope are plotted in the same format as in Figure El here relative to 
the position at t = 123rA. 



fast ascent (see next section). Despite the very strong distortion, the flux rope returned close 
to the original position, executing quickly decaying relaxation oscillations. This is shown in 
Figure [9] for the perturbation applied at t = 123rA. Nearly identical behavior results for the 
perturbation applied at the end of the relaxation run. These tests support the interpretation 
that the model with axial flux $axi = 5 x 10^^ Mx relaxes to a stable equilibrium. 

The model can be driven to eruption by further flux cancelation. This was tested by 
prescribing horizontal flows converging toward the PIL in the strong-flux section of the PIL 
a nd allowing for diffusio n of the fleld in the vic inity of the PIL, analogous to the diflFusion 



m lAulanier et al.l ( 120101 ) and lAmari et al.l ( 120111 ). Cancelation of about 10% of the unsigned 
flux in the region caused the eruption of the flux rope (but more extended experimenting 
might flnd eruption at an earlier stage). DiflFerent from observation, the flux rope took a 
nearly vertical path. This and the fact that an appropriate modeling of pre-eruption driving 
by flux cancelation should start from a magnetogram taken at an earlier time, led us to 
reserve a more detailed investigation of such simulations for a future study. 
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6. Eruption of the Unstable Configuration 

Similar to the result of the MFR in Paper I, the rope with axial flux $axi = 6 x 10^^ Mx 
does not relax to an equilibrium but erupts. The evolution starts with a rapid displacement 
of the rope within the flrst two Alfven times, comprising the lifting of the middle part and the 
southeastward expansion of the eastern elbow which are very similar to the initial behavior 
in the relaxation runs of the previous section. The interaction with the flux rooted in the 
minor positive polarity is also similar but does not lead to a split rope (which occurs at 
t > 5ta in the relaxation runs) because of the commencing eruption. 

A seamless transition to an initially slower but accelerating rise follows, leading to the 
eruption of the whole flux rope (see Figure [10] and the accompanying animation) . The main 
direction of this motion is upward and southward, initially (at t = 2ta) at an inclination 
of 40°-45° to the ver tical, which is very close to the observed direction. Interaction with a 



nearby co ronal hole (|GoDalswamv et al.l l2009l ) and asymmetry of the photospheric flux dis- 
tribution ( jPanasenco et al.ll2012l ) have been proposed to cause deviations from radial ascent. 



The former is not included in our model, while the latter eflFect is present, but presumably 
not very strong: the main positive polarity in the active region contains only ^ 10% more 
flux than the main negative polarity. However, the positive polarity is more compact. An 
additional effect is suggested by the shape of loop sets 1 and 3: a southward-directed sling- 
shot action by their common flux after the flux rope loses equilibrium. The fact that the 
dimmings of the eruption begin to develop near their end points (Figure [11(d)) indicates 
that these loop sets possess common flux which plays a role in the eruption. DiflFerent from 
observation, the rise soon begins to gradually turn more radial (at t = 6ta the inclination 
is ^ 30°). This may be due to the missing action of the polar coronal hole and due to the 
slingshot eflFect being weaker than in reality, since the time for loop set 3 to develop fully is 
not available in the unstable model. The latter aspect may be improved by applying longer 
magnetofrictional relaxation in the preparation of this model. 

The outer layers of the expanding flux rope begin to experience compression at the side 
boundary of the computation box already at ^ 6rA, while the core hits the boundary with 
a delay of (l-2)rA. The evolution up to this point shows a number of similarities to the 
observed behavior (see also below). Subsequently, the upward deflection of the rope is of 
course unrealistic, but the continuing strong expansion of the minor flux rope radius and the 
continuing reconnection under the rope, discussed below, are likely qualitatively consistent 
with the modeled solar event — until the top boundary is encountered. 

The eruption of the rope lifts the overlying flux. This flux is initially sheared (Figure [TOl) 
but becomes more antiparallel in the volume under the rope as it is lifted. The corresponding 
induced current points primarily horizontally and is associated with a Lorentz force pointing 
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toward the essentially vertical layer under the flux rope where the vertical fleld component 
of the ambient flux changes sign. The induced current thus pinches into the vertical current 
sheet, or flare current sheet, known from the standard flare model. One can understand this 
pinching also as the 3D generalization of the well-known instability of a magnetic X point 
in two dimensions. The generalized X-type structure is known as a hyperbolic flux tube 
(HFT). Its pinch ing into a current shee t following an appropr iate perturbation has been 
demonstrated by iTitov et al.l ( 120031 ) and iGalsgaard et al.l ( 120031 ) . The initial conflguration 
contains such a structure (see Figure [3] and Paper I) . Note that the pinching is solely driven 
by the Lorentz force in our zero-beta simulation (the pressure gradient would additionally 
contribute if /3 > 0). 

Figure [TT] displays the vertical current sheet that forms under the rising flux rope. Since 
it is a true current sheet (with exponentiall y rising current de nsity and correspondingly de- 
creasing current sheet width in ideal MHD; iTitov et al.ll2003l ). the pinching process reaches 
saturation in a standard one-fluid MHD description only if suflicient diffusion is provided. 
Otherwise, it proceeds toward the limit of the employed numerical scheme, resulting in un- 
physical fllamentation of the current sheet (w hich develops on the grid scale and is not re lated 
to the plasmoid instability of current sheets; iLeboeuf et al.lll982l : iLoureiro et al.l 120071 ). To 
verify the eruption of the third model under uniform numerical settings, and as a reference, 
we have flrst followed the latter evolution in a run without magnetic fleld smoothing, a^ = 0. 
Magnetic reconnection develops in the vertical current sheet and forms coherent reconnection 
outflows in spite of the ensuing fllamentation, so the main ingredients of an eruption — loss 
of equilibrium and magnetic reconnection — are present in spite of the numerical artefacts. 
Repeating this run with magnetic fleld smoothing, chosen to be uniform in the flrst instance, 
we flnd that the flux rope rise velocity is not strongly influenced as long as a^ < 10~^-^; it 
then stays within 20% of its nominal value for cr^ = 0. Beyond that level, the rise velocity 
falls off quickly, and for as > 10~^-^ all velocities diffuse away after the initial upward dis- 
placement of the rope. It is clear that when the magnetic diffusion as is large, the currents 
and Lorentz forces in the body of the flux rope are weaker, artiflcially stabilizing the solution. 
The fllamentation of the current sheet decreases with increasing as but disappears only for 
(^B ^ 10~^. Therefore, the appropriate numerical setting for the present model consists in 
applying the smoothing of the magnetic fleld only in the volume under the flux rope using 
a level of a^ = 10~^, with a gradual transition to a^ = at the bottom and at the sides. 
This avoids the fllamentation of the vertical current sheet and yields a rise velocity of the 
flux rope apex very close to (slightly higher than) the value obtained in the absence of the 
smoothing. Figures [T0l - [T6l show the results of this run. 

In addition to the pinching of the HFT under the rising flux rope, part of the current 
layer at the periphery of the rope steepens, while the main part of the current layer weakens 
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(since the total current through the rope must decrease to power the eruption). The bottom 
part of the current layer at the side of the flux rope facing the stronger ambient fleld north- 
eastward of the PIL builds up current densities comparable to the ones in the pinched HFT 
underneath; thus, the vertical current sheet extends upwards asymmetrically along this side 
of the flux rope. 

The flows triggered by the eruption are shown in the second panel of Figure dH One 
can see the rise and expansion velocity in the whole cross section of the flux rope, as well as 
in the ambient volume, which is threaded by fleld lines lifted by the rope. The southward 
inclination of the eruption is apparent. Reconnection flows develop at the vertical current 
sheet essentially simultaneous with the rise of the rope. The flrst indication of the downward 
reconnection outflow becomes visible at t ?^ Ita, i.e., already in the initial upward displace- 
ment to the rope's equilibrium height. As the eruption progresses, the reconnection flows 
evolve joint with the vertical current sheet; thus, they become aligned with the northeast- 
ward periphery of the rope in the mai n phase of the eruption shown in the flgure. This may 



yield a contribution to the roll effect (jPanasenco et al.ll201ll ). which is weakly indicated by 



animated fleld line plots from this simulation, as well as in the EUVI-A data shown below 
in Figure [151 Finally, we note that the upward reconnection outflow reaches slightly higher 
velocities than the flux rope and its immediate surroundings, so that the new flux carried 
into the rope is likely to have a direct positive eflFect on the acceleration in the main body 
of the rope. 

It is worth noting that the stable model with $axi = 5 x 10^° Mx also possesses an HFT 
(see Figure [3](c) above and Figure 7(d) in Paper I). This HFT pinches as well into a short 
vertical current sheet during the initial upward displacement, and reconnection commences 
as early as for the unstable model; see the direct comparison of the models at t = 2tx 
(i.e., after the initial displacement) in Figure [T2l However, an eruption does nevertheless 
not occur. This indicates that reconnection in the vertical current sheet of the considered 
coronal NLFFF models is not able to drive an eruption by itself, but rather that the driving 
by an unstable flux rope is required. A downward return flow is induced in the stable case 
by the rapid initial rise of the flux rope (Figure [T2l(a)): this would not occur in a quasi-static 
evolution. This flow is seen to join the reconnection inflow, so it does not appear to work 
against the further development of reconnection in this model. 

Figure [13] compares flare loops and reconnected fleld lines in the simulation. The selected 
observation times correspond to t = 2.8rA and t = 7.9rA in the simulation if the scaling of 
Section 16.11 is adopted. A similar comparison with fleld lines of the corresponding MFR 
model is given in Figure 12 of Paper I. The two sets of green and red fleld lines are traced 
from start points at z = in the positive polarity which lie slightly inside the current layer 



22 



that extends from the bottom tip of the vertical current sheet at 2.8 and 7.9rA, so that they 
are to be compared with the flare loops in panels (c) and (d), respectively. (In the flgure, 
both sets are traced in the fleld at 7.9ta] the lower set is visually nearly indistinguishable 
from the fleld lines traced in the fleld at 2.8rA.) The trend of progressively decreasing shear 
is very clear in the simulation data as well, see the vertical view in panel (a). The perspective 
view at an inclination of 25° from vertical in panel (b) is chosen to correspond to the latitude 
of ^ 25° of the active region (the weak tilt corresponding to the slightly eastward longitude 
is not taken into account). Here the fleld line shapes appear basically similar to the shapes 
of the flare loops. The loops possess slightly higher shear. This is not unexpected, since the 
use of the potential fleld in constructing the NLFFF models is likely to remove some of the 
shear that has accumulated in the ambient fleld in the earlier evolution of the active region. 
Part of this shear is recovered in the MFR phase, and the remaining diflFerence appears to 
be minor for the considered region. 

A more signiflcant and also expected diflFerence is revealed by the distance of the flare 
ribbons to the PIL. Compared to the end points of the newly reconnected fleld lines at 7.9rA, 
the separation of the flare ribbons at 03:23 UT is larger by a factor ?^ 1.5, which signifles 
a corresponding diflFerence in the reconnection rate. Given that no attempt was made to 
correctly model the reconnection rate (beyond the choice of the numerical diffusion such 
that the current sheet does not develop numerical artefacts), the difference is surprisingly 
small. We expect that it can be removed by including resistivity. 

The rise proflle of the eruption is plotted in Figure [HI Here the initial, mainly upward- 
directed displacement leads to a doubling of the initial apex height (from ho = 0.047 to 
h{2) = 0.098) within a time frame similar to the relaxation runs. The flux rope has roughly 
semicircular shape at this time. As soon as the initial velocity enhancement has largely 
decayed, one can see a more gradually developing rise with a signiflcant horizontal component 
(obvious from the increasing difference between the solid and dashed curves in the flgure). 
The logarithmic display reveals that the more gradual rise is approximately exponential 
(but it can also be approximated by a power law with an exponent near 3). The different 
functional form and direction of the rise after t = 2ta suggest that the two phases have 
different physical origins. We interpret the initial upward displacement as the motion to the 
equilibrium position of the flux rope which it had not yet reached after the 30,000 MFR 
iterations that deflne the initial condition for the MHD run. The subsequent exponential-to- 
power law rise is charac teristic of an instability launched by a small-to-moderate perturbation 



( jSchrijver et al.ll2008bl ) from the obviously unstable equilibrium position. 



The fluid element at the magnetic axis of the flux rope experiences the reflection at the 
side boundary of the box after t ^ 8ta- 8.7 (39) percent of the initial total (free) magnetic 
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energy are released at this stage. The subsequent part of the rise profile is not related to 
reality and included only for completeness. The additional upward acceleration in the process 
of the reflection results from the induction of currents as the expanding flux is compressed 
at the side boundary. 
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Fig. 10. — Flux rope eruption in the model with $axi = 6 x 10^° Mx. The format is similar 
to Figs, m and [71 The whole box is shown in the left two columns. The eruption starts 
strongly inclined, as observed. Subsequently, the erupting flux is deflected at the closed side 
boundary. Originally overlying flux (green fleld lines) is strongly reconnected into the rope 
by flare reconnection in the vertical current sheet that forms under the rope (see Figure fTTI) . 
(A color version and an animation of this flgure are available in the online journal.) 
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Fig. 11. — Vertical current sheet in the unstable configuration at t = 6.6ta. {Left): Isosurface 
of current density at 0.1 max(|j |). (Right): Force- free parameter a{s^ z) and in-plane velocity 
vectors in the vertical cut plane shown in Fig. [3l 
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Fig. 12. — Reconnection fiows of the models with (a) $axi = 5 x 10^° Mx and (b) $a 
6 X 10^° Mx at t = 2rA. 
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Fig. 13. — Reconnected field lines in the unstable configuration in (a) vertical view and (b) 
perspective view inclined by 25°. (c-d) Flare loops of the event at the times corresponding 
to, respectively, t = 2.8ta and t = 7.9ta in the simulation. Contours of the vertical field 
component are included in (c). (A color version of this figure is available in the online 
journal.) 
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Fig. 14. — Rise profile of the flux rope apex in the unstable conflguration, shown in the 
left panels in a format similar to Fig. [6l The right panels show the same information on a 
logarithmic scale. The section of the rise after reflection at the side boundary is plotted in 
gray. 
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6.1. Scaling the Simulation to the Observations 

From the rise profile in Figure [T^ it is obvious that the simulation can capture properties 
of the eruption on 2010 April 8 in the early, approximately exponential phase, 2 < t/rx < 8. 
The quantitative comparison between this part of the rise profile and EUVI-A 195 A data 
allows us to scale the simulation to the observed event. Since the erupting flux propagates 
mainly southwards (as seen in the AIA images), STEREO Ahead observes its rise nearly 
exactly from the side, presumably with relatively small projection effects. Note that only 
the scaling of the time unit ta, or equivalently of the velocity unit Vao, is to be determined, 
since the length and field strength units are given by the magnetogram and the density is 
fixed by field strength and Alfven velocity. 

An animation of the EUVI-A 195 A images is included with Figure [151 ^^d panel (a) 
shows a representative image. No trace of the rising filament material can be seen, but 
several overlying loops stand out relatively clearly. Their rise is quantified in the stack plot 
(panel (b)) taken at an artificial slit which is aligned with the average direction of ascent 
(see panel (a)). The three clearest traces are marked in color. 

Next, we select field lines which correspond to the overlying loops. It turns out that 
no exact match is possible, since the usable height range in the simulation falls between the 
lower two traces in the stack plot. This is seen in Figure [T5l(c-d). where a ray is indicated 
which corresponds to the artificial slit in the EUVI-A image. The panels display a view in 
—X direction, similar to the view of STEREO Ahead. The ray starts at the apex of the flux 
rope's magnetic axis at t = and lies in the y-z plane, inclined to the vertical by 45°. After 
the initial upward displacement, the flux rope radius has grown such that the innermost 
overlying flux lies at a distance of about 70 Mm along the ray, larger than the pre-eruption 
distance of ^ 40 Mm of the green trace in the stack plot. The pre-eruption distance of the 
red trace is about 175 Mm, already relatively close to the side boundary of the box. Flux 
in this height range can be followed until the compression zone near the side boundary is 
encountered, roughly at 250 Mm along the ray, but this is of course stronger influenced by 
the side boundary than the flux immediately surrounding the rope, and it can utilize only a 
small part of the red trace in the stack plot. Nevertheless, using a field line in the innermost 
overlying flux permits an acceptable scaling, since the expansion of the overlying flux is 
rather uniform: the three traces run approximately parallel to each other, i.e., with similar 
velocity profiles. 

We thus select a fluid element on the ray which lies slightly outside the flux rope at 
t ^ 1.75rA, the beginning of nearly exponential rise. The trajectory of this fluid element 
is computed through the simulation. At later times, a field line is traced from the new 
position and projected onto the y-z plane, where its crossing with the ray yields a new 



29 



projected distance along the ray. Finite differencing yields a projected velocity along the 
main direction of propagation for both observation and simulation data. 

The EUVI-A images (of 2.5 min cadence) place the onset time of the expansion in 
a relatively narrow interval, 02:30:30-02:33:00 UT, of which we associate a time near the 
midpoint, 02:32 UT, with t = 1.75ta in the simulation. The value Vao ^ 1400 kms~^ that 
yields a good match of both the projected distance and velocity data is then easily found by 
trial and error; see the diamonds in Figure [T5l(e-f). which cover the interval t = (1.73-5.93)rA 
in the simulation. A change of ±100 kms~^ from the estimated Vao already produces a much 
poorer match. 

We note that a trend of decreasing acceleration becomes visible in the projected velocity 
data of the overlying flux in the simulation already in the middle of the considered interval. 
This is not seen in the observation data and indicates how early the overlying flux is influenced 
by the side boundary. Nevertheless, the match between simulation and observation appears 
acceptable up to the second-to-last velocity point in the interval. 

An ambiguity in the above procedure originates from the insuflicient knowledge about 
the position of the overlying loops relative to the plane of sky for STEREO Ahead. They 
cannot be unambiguously identifled in the corresponding AIA 193 A images. We know 
from AIA that the main direction of expansion was nearly southward, i.e., nearly in the 
plane of sky for STEREO ALead, but that the initial expansion had a strong southeastward 
component. To roughly estimate how this ambiguity might influence the scaling, we repeat 
the procedure, positioning the fluid element at t ?^ 1.75rA on rays directed ±30° off the 
meridional plane (but also starting at the flux rope apex). For the southeastward pointing 
ray, the same Vao is obtained from a match of nearly the same quality (see the plus signs in 
Figure [T5l(e-f)). For the southwestward pointing ray, a higher Alfven velocity is indicated (by 
^ 50%), but it is not possible to match both distances and velocities to the stack plot data 
at a reasonable quality. The AIA data do not show any distinctive structure propagating 
in this direction in the time interval used in the scaling. Thus, we consider it very unlikely 
that the structures seen by EUVI-A have expanded in this direction. On the other hand, the 
consistency of the Vao estimates, based on the southward and southeastward directed rays, 
makes them rather reliable. For completeness we note that the considered fluid elements 
remain relatively close to their respective ray in the course of the eruption. 

The estimated value of the peak initial Alfven velocity Vao refers to the point of peak 
field strength in the box, max(So) = 955 Gauss, which lies in the magnetogram plane. Our 
model for the initial density implies VA(t = 0) oc Bq^ , so that the initial field strength of 
^ 70 Gauss in the area of the fiux rope apex yields an Alfven velocity VA{t = 0) ^ 730 kms~^ 
and a particle density N{t = 0) ^ 4.4 x 10^^ cm~^ in the upper part of the fiux rope. The 
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density must be taken as a characteristic average density in the filament channel, since our 
density model disregards any filamentary fine structure of the plasma. 

One has to consider this velocity as the lower end point of a plausible range (and the 
density, correspondingly, as an upper end point). This value follows from a model with 
a certain amount of prescribed axial flux, $axi = 6 x 10^° Mx, which lies near the point 
of marginal stability, but at a somewhat arbitrary distance set by the interval of 10^^ Mx 
between the $axi values in consecutive models. An unstable model with less axial flux would 
yield a slower rise in the simulation, consequently, a higher estimate for Vao- Indeed, the 
leading-edge CME velocity of ^ 520 km s~^ observed by C0R2 is at 70% of the Alfven 
velocity in the source. Allowing for a factor 2 lower velocity of the CME core, this becomes 
^VA(t = 0)/3, still a rather high value for a moderate CME that gains its acceleration in 
a large height range of sever al Re?). Such eruption s tend to reach lower terminal flux rope 
velocities of < 0.2VA(t = 0) (JTorok fc Klieml 120071 ) . Thus, a range for the Alfven velocity 
in the flux rope from ^ 730 kms~^ up to about a factor 2 higher appears reasonable. The 
corresponding density range is ^ (1^4) x lO"*^^ cm~^. 



These estimates lie well within acceptable ranges. ICaryl (1200 ll ) finds /3 ^ 0.02-0.2 at a 
height of 0.1 Rq (see his Figure 3), which translates to Va = IC^^'^I^ - (650-2000) km s"^ 
for a plasma temperature of T = 2.5 MK (i.e., a sound speed of C^ = 145 km s~^). The 
lower part of the ran ^e is appropria t e for a considerably dispersed, moderate active re- 
gion like AR 11060. iLabrosse et al.l ( 120101 . their Section 3.3.2) quote the wide range of 
A^ ~ (10^-10-^-^) cm~^ for densities of cool (Ha-emitting) to moderately hot (EUV-emitting) 
plasmas in quiescent and active-region prominences. The lower-to-middle part of the range 
may be representative for the average density in the erupting flux in the considerably dis- 
persed active region. 
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Fig. 15. — Scaling of the simulation to the observation data, (a): A representative EUVI-A 
195 A radial-filter image showing the structure of the overlying flux. The indicated artificial 
slit is aligned with the main expansion direction of 45° from vertical and used to generate 
the stack plot in panel (6). (c-d): Field line plots during the nearly exponential rise phase, 
viewed in —x direction. The field line used for the scaling is marked in blue, (e-/): Scaling 
of simulation data (blue diamonds and plus signs) to the rise profiles marked in the stack 
plot. (A color version of this figure and an animation of the EUVI-A 195 A image sequence 
are available in the online journal.) 
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6.2. Onset Condition 



By reaching good agreement with key features of the eruptive event on 2010 April 8 
hke the initial height-time profile and the rise direction of the erupting fiux, our simulations 
substantiate the conclusion in Paper I that the fiux rope insertion method allows to model 
the NLFFF in AR 11060 around the time of the eruption. The condition for the loss of 
equilibrium, formulated as a condition of fiux imbalance in terms of the ratio between the 
rope's axial fiux and the unsigned fiux in the source region of the eruption, is also confirmed. 
Given the total unsigned fiux in AR 11060 of F^ = 3.7 x 10^^ Mx (Paper I), the range of 
^axi = (5-6) X 10^° Mx for the marginal stabilit y point yields a flux ratio ^r,^-J(F,j2) = 



(27-32)%, well within the range found previously (JBobra et al.ll2008l : ISavcheva et al.ll2012bl ). 
A similar flux ratio is reached at the onset of eruption in the simulation that applied flux 
cancelation to the model with $axi = 5 x 10^^ Mx. 

It is of principal interest to compare this condition with the condition obtained by 
describing the loss of equilibrium as a plasma instability, i.e., the kink instability of a current 
channel. The relevant mode of the kink instability here is the lateral kink, known as the 
torus instability in the case of an arched current channel, since the pre-eruptive NLFFF 
possesses insufficient twist for the excitation of the helical kink mode. Using the length and 
fiux values of the two fiux rope models next to the marginally stable point, L ^ 270 Mm, 
$axi = (5-6) X 10^° Mx, and Fpoi = 10^° Mx cm"^ the average twist angle is LB^/{rB^) ^ 
7rLFpoi/$axi = (0.53-0.45)7r, far smaller than the threshold of the helical kink. 

The threshold of the torus instability is given in terms of the decay index of the external 
poloidal field, Sgp, at the position of the current channel, n = — rflnSgp/rflnz > ncr. The 
canonical values of the critical decay index are n^ = 3/2 fo r a toroidal current channe l 
( BatemanI 1 19781 ) and rirr = 1 for a straight current channel (Ivan Tend fc KuperusI Il978l ). 
Kliem fc Torokl (120061 ) showed that the value for a toroidal current channel rises to ricr = 2 if 
reconnection under the fiux rope does not commence im n iediat ely, and they added a relatively 
small negative correction term. iDemoulin fc Aulanierl (J2010l ) generalized the consideration 
to arbitrarily shaped paths of the current channel, finding a range ricr ^ 1.1-1.3. Most 



nume rical studies, inc 



1.75 (Torok fc Kliem 



udin g the only parametric study to d ate, find ricr in the range 1.5 



20071: lAulanier et all l2010l : iFanI l2010l ). but it should be noted that 



the simulation in iFan fc GibsonI (l2007l) yields a value of 1.9. The range 1.5-1.75 is also 
supported by observational studies (JLiul l2008l : IGuo et al.l I2OIOI ) . The critical decay index 
must depend upon the photospheric line tying (which varies with the shape of the current 
channel) and it should increase for increasing strength of the external toroidal (shear) field. 
These dependencies have not yet been investigated. However, we can expect these effects 
to be relatively minor in the present case, since the line tying is minimized for semicircular 
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flux rope shape, which the unstable flux rope reaches at the onset of the eruption, and since 
the ambient potential fleld has only little shear. In the comparison below, we adopt the 
currently best supported range of the torus instability threshold, n^ ^ 1.5-1.75. 

In practice it is often impossible or rather laborious to isolate the external poloidal fleld 
Sep. This involves the subtraction of the poloidal fleld created by the current channel in the 
corona, which requires the exact knowledge of the current channel, i.e. of the NLFFF. Here 
we have this knowledge from the relaxed conflguration with $axi = 5 x 10^^ Mx and from 
the initial conflguration with $axi = 6 x 10^° Mx. However, given the effort required for this 
analysis and given the uncertainty in the knowledge of the critical decay index, we choose to 
simplify the consideration by instead using the decay index of the horizontal component of 
the (total) potential fleld, h = —dlnB^ot^xy/d^^^- We have verifled that the potential fleld 
passes at nearly right angles over the PIL in the height range of the flux rope at t = 0, so that 
it is a reasonable approximation of the external poloidal fleld. The direction is much closer 
to perpendicular than that of the immediately overlying flux in the initial conflguration of 
the MHD simulation (green fleld lines in Figure [TOl at t = 0), which is considerably influenced 
by the presence of the flux rope in the course of the MFR. We have computed h{z) at several 
positions along the PIL and found it to be nearly uniform in the section of the PIL between 
the main flux concentrations where the middle part of the inserted flux rope is located. 
Figure [16] shows the proflle at the vertical line through the flux rope apex. The height range 
corresponding to the adopted range of ricr is z = 0.063-0.094. 

The bottom panel of Figure [16] summarizes the height-time proflles of the three models 
studied in this paper. The apex of two stable ones stays in the range h ^ 0.03-0.055, i.e. 
below the torus-unstable height range. The inferred equilibrium apex height of the unstable 
model, h{t = 2) = 0.098, falls rather clearly in the torus-unstable height range. Thus, for 
the considered AR 11060, the stability boundary given in terms of the axial flux in the rope, 
$axi = (5-6) X 10^° Mx, corresponds very well to the most likely range for the threshold of 
the torus instability, ricr ^ 1.5-1.75. 
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Fig. 16. — (a) Decay index h = —dlnBpot,xy/dlnh of the horizontal components of the 
potential field along a vertical line passing through the apex point of the fiux rope's magnetic 
axis at t = 0. (b) Height-time profiles of the fiux rope apex for the three models studied in 
this paper, compared to the onset condition of the torus instability. (A color version of this 
figure is available in the online journal.) 
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7. Summary and Conclusions 



This paper reports MHD simulation s that study NL FFF models of the eruptive AR 11060 



The models were constructed in Paper I (jSu et al.ll201ll ) using the flux rope insertion method 



and magnetofrictional relaxation. We selected three that appeared to best represent the pre- 
emptive, nearly marginally stable, and eruptive states of the active region, differing only in 
the axial flux of the inserted flux rope, which encompasses the range $axi = (4-6) x 10^^ Mx. 
Our simulations conflrm key results of Paper I: 

1) There is a limiting value of the axial flux in the rope for the existence of stable NLFFF 

equilibria, which lies in the selected range. The equilibrium height of the flux rope 
increases with axial flux. 

2) Stable models relax deeply in the MHD simulation, reaching numerical equilibria that 

retain a flux rope in the considered range of $axi- Field lines in the flux rope and its 
immediate surrounding correspond well in shape to the observed coronal loops. 

3) The flux rope in the unstable model experiences a full eruption. 

Thus, although none of the observed structures in AR 11060 shows regular multiple crossings, 
the deflnite signature of a flux rope, our simulations substantiate the conclusions in Paper I 
that a flux rope existed in the active region prior to the eruption and that the rope's loss of 
equilibrium caused the eruption. Additionally, we obtain the following results: 

4) The model with $axi = 5 x 10^° Mx, found to be nearly marginally stable in Paper I, is 

also found to be closest to marginal stability — on the stable side. It does not erupt even 
though the hyperbolic flux tube at its underside pinches into a short vertical current 
sheet and a transient phase of reconnect ion in this current sheet is triggered, lasting 
for a couple of Alfven times. The model can be driven to eruption by photospheric 
flux cancelation, and the ratio of its axial and unsigned flux at that point is similar to 
the flux ratio of the unstable model. 

5) The erupting flux rope shows an accelerated rise /i(t), which can be approximated by an 

exponential and also b y a power-law with an i ndex near 3. Such a rise is characteristic of 



a flux rope instability (jSchrijver et al.ll2008bl ). The rise velocity reaches a considerable 



fraction of the initial Alfven velocity in the flux rope, about 20% (still rising when 
the boundary of the simulation box is approached). This indicates that the model lies 
clearly in the unstable domain of parameter space. 
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6) The simulated eruption agrees very well with several characteristics of the modeled event, 

including (i) the initially very inclined direction of ascent, by 45° from vertical toward 
the equator; (n) the accelerated rise of the overlying flux; and (m) the close temporal 
association between the acceleration of the flux rope (the CME component of the 
eruption) and the development of reconnection in the vertical current sheet underneath 
(the flare component). Additionally, initially high shear and the trend of progressively 
decreasing shear in the reconnected flux under the rope, as displayed by the flare 
loops, are clearly reproduced. However, there are also differences, (i) The simulated 
rise begins to turn more radial already low in the corona, whereas the modeled event 
showed an increasing inclination in the low corona, turning more radial only at heights 
of several solar radii, (ii) The footpoint locations of newly reconnected fleld lines in 
the simulation move away from the PIL at a somewhat lower speed than the observed 
flare ribbons, indicating a moderately lower reconnection rate. Obvious reasons and 
straightforward options for improvement of the modeling were identifled for both. 

7) Oblique eruption paths can be caused by the speciflcs of the magnetic structure in the 

source region, in addition to the interaction with a coron al hole and the asymme- 



try of the photospher ic flux distribution suggested earlier (JGopalswamy et al.l l2009l : 



Panasenco et al.l 1201 2l ) . We suggest that the slingshot action of horizontally curved 
flux, which was or quickly became part of the erupting flux, contributed to the incli- 
nation of the rise path in the modeled event. 

8) The ability to follow the eruption allows us to scale the simulation to the observation 

data, thus estimating the Alfven velocity in the source region of the eruption. The 
range Va ^ (730-1500) kms~^ in the volume of the flux rope is indicated. Since the 
fleld strength in the rope is given by the NLFFF model, B ^ 70 Gauss, the estimated 
Alfven velocity implies an average flux rope density in the range A^ ^ (1~4) x 10^° cm~^. 
Both values lie within acceptable ranges, especially the Alfven velocity, since AR 11060 
was already considerably dispersed. 

9) By its sharper discrimination between stability and instability compared to the mag- 

netofrictional relaxation technique, the MHD modeling yields a narrower range for the 
threshold value of the considered control parameter. We flnd that the value lies in the 
range $axi = (5-6) x 10^° Mx (compared to (4-6) x 10^° Mx in Paper I). This corre- 
sponds to a range of the flux ratio $axi/(i^u/2) = 0.27-0.32, where F^ = 3.7 x lO^^ Mx 
is the total unsigned flux in the active region. 

10) This range of the flux ratio corresponds very well to the threshold of the torus instability 

in the considered active region. Using the horizontal component of the potential fleld as 
a proxy for the external poloidal fleld component, the decay index n = — rflnSgp/rflnz 
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takes values n ^ 1.3^1.8 in the range of equilibrium flux rope heights found for the 
range of axial flux $axi = (5-6) x 10^° Mx. The critical decay index for onset of the 
torus instability is estimated to lie in the range n^ ^ 1.5-1.75. This establishes a 
connection between these independently developed criteria for the loss of equilibrium 
of a flux rope, which is plausible from the following consideration. Higher axial flux 
in the rope corresponds to greater equilibrium height (as a consequence of the implied 
higher flux rope current), which, in turn, typically corresponds to higher decay index, 
i.e., to approaching the torus-unstable height range where n > ricr- 

11) A considerable amount of changes was found while the stable models relaxed to a 
numerical equilibrium in the MHD simulation, although they had already gone through 
a magnetofrictional relaxation phase. The magnetofrictional relaxation had a pre-set 
number of iterations, thus it had not progressed to the deepest level possible, especially 
not for the nearly marginally stable model, but the model with $axi = 4 x 10^^ Mx 
was quite well relaxed. Even changes of the topology (a transient split of the flux 
rope) occurred in the MHD relaxation. Only part of the different level of dynamics 
appears to be due to the numerical differences between the employed magnetofrictional 
and MHD simulation codes. A systematic comparison of magnetofrictional vs. MHD 
relaxation might thus be warranted. 



Based on various developments and applications of t he flux rope insertion method (e.g. 



van Ballegooijenll2004l : iBobra et all 12008 : 



Su et al.ll201ll ). this research represents a further 



step toward data-constrained numerical modeling of solar eruptions. Future work along this 
line will consider coronal NLFFF const ructed from a t i me ser ies of magnetograms, as in 
Savcheva fc van BallegooijenI ( 120091 ) and ISavcheva et al.l (l2012bl ) , and eventually proceed to 
driving coronal NLFFF through time-dependent photospheric boundary data obtained from 
observations. Technical improvements toward using larger domains and open boundary 
conditions will permit following eruptions much further than possible here; this will set 
tighter constrains as well. 
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